﻿#region Sharp3DMath, Copyright(C) 2004 Eran Kampf, Licensed under LGPL.
//	Sharp3DMath math library
//	Copyright (C) 2003-2004  
//	Eran Kampf
//	eran@ekampf.com
//	http://www.ekampf.com
//
//	This library is free software; you can redistribute it and/or
//	modify it under the terms of the GNU Lesser General Public
//	License as published by the Free Software Foundation; either
//	version 2.1 of the License, or (at your option) any later version.
//
//	This library is distributed in the hope that it will be useful,
//	but WITHOUT ANY WARRANTY; without even the implied warranty of
//	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
//	Lesser General Public License for more details.
//
//	You should have received a copy of the GNU Lesser General Public
//	License along with this library; if not, write to the Free Software
//	Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
#endregion
#region Using directives

using System;
using System.Drawing;
using System.Drawing.Drawing2D;
using System.Collections.Generic;
using System.Text;
using Sharp3DMath.Strokes;
using Sharp3DMath.Geometry2D;



#endregion

namespace Sharp3DMath.Core
{
	/// <summary>
	/// Provides standard mathematical functions for the library types.
	/// </summary>
	public static partial class MathEx
	{
		#region Delegates
		public delegate T VoidFunction<T>();
		public delegate T UnaryFunction<T>(T param);
		public delegate T BinaryFunction<T>(T param1, T param2);
		public delegate T TenaryFunction<T>(T param1, T param2, T param3);
		#endregion

		#region Constants
		/// <summary>
		/// The value of PI.
		/// </summary>
		public const double PI = System.Math.PI;
		/// <summary>
		/// The value of (2 * PI).
		/// </summary>
		public const double TwoPI = 2 * System.Math.PI;
		/// <summary>
		/// The value of (PI*PI).
		/// </summary>
		public const double SquaredPI = System.Math.PI * System.Math.PI;
		/// <summary>
		/// The value of PI/2.
		/// </summary>
		public const double HalfPI = System.Math.PI / 2.0;

		/// <summary>
		/// Epsilon, a fairly small value for a single precision floating point
		/// </summary>
		public const float EpsilonF = 4.76837158203125E-7f;
		/// <summary>
		/// Epsilon, a fairly small value for a double precision floating point
		/// </summary>
		public const double EpsilonD = 8.8817841970012523233891E-16;
		#endregion

		#region Function Definitions

		#region Trigonometric Functions
		public static readonly UnaryFunction<double> DoubleSinFunction = new UnaryFunction<double>(MathEx.Sin);
		public static readonly UnaryFunction<double> DoubleCosFunction = new UnaryFunction<double>(MathEx.Cos);
		public static readonly UnaryFunction<double> DoubleTanFunction = new UnaryFunction<double>(MathEx.Tan);
		public static readonly UnaryFunction<double> DoubleCotFunction = new UnaryFunction<double>(MathEx.Cot);
		public static readonly UnaryFunction<double> DoubleSecFunction = new UnaryFunction<double>(MathEx.Sec);
		public static readonly UnaryFunction<double> DoubleCscFunction = new UnaryFunction<double>(MathEx.Csc);

		public static readonly UnaryFunction<ComplexD> ComplexDSinFunction = new UnaryFunction<ComplexD>(MathEx.Sin);
		public static readonly UnaryFunction<ComplexD> ComplexDCosFunction = new UnaryFunction<ComplexD>(MathEx.Cos);
		public static readonly UnaryFunction<ComplexD> ComplexDTanFunction = new UnaryFunction<ComplexD>(MathEx.Tan);
		public static readonly UnaryFunction<ComplexD> ComplexDCotFunction = new UnaryFunction<ComplexD>(MathEx.Cot);
		public static readonly UnaryFunction<ComplexD> ComplexDSecFunction = new UnaryFunction<ComplexD>(MathEx.Sec);
		public static readonly UnaryFunction<ComplexD> ComplexDCscFunction = new UnaryFunction<ComplexD>(MathEx.Csc);
		#endregion

		#region Trigonometric Arcus Functions
		public static readonly UnaryFunction<double> DoubleAsinFunction = new UnaryFunction<double>(MathEx.Asin);
		public static readonly UnaryFunction<double> DoubleAcosFunction = new UnaryFunction<double>(MathEx.Acos);
		public static readonly UnaryFunction<double> DoubleAtanFunction = new UnaryFunction<double>(MathEx.Atan);
		public static readonly UnaryFunction<double> DoubleAcotFunction = new UnaryFunction<double>(MathEx.Acot);
		public static readonly UnaryFunction<double> DoubleAsecFunction = new UnaryFunction<double>(MathEx.Asec);
		public static readonly UnaryFunction<double> DoubleAcscFunction = new UnaryFunction<double>(MathEx.Acsc);

		public static readonly UnaryFunction<ComplexD> ComplexDAsinFunction = new UnaryFunction<ComplexD>(MathEx.Asin);
		public static readonly UnaryFunction<ComplexD> ComplexDAcosFunction = new UnaryFunction<ComplexD>(MathEx.Acos);
		public static readonly UnaryFunction<ComplexD> ComplexDAtanFunction = new UnaryFunction<ComplexD>(MathEx.Atan);
		public static readonly UnaryFunction<ComplexD> ComplexDAcotFunction = new UnaryFunction<ComplexD>(MathEx.Acot);
		public static readonly UnaryFunction<ComplexD> ComplexDAsecFunction = new UnaryFunction<ComplexD>(MathEx.Asec);
		public static readonly UnaryFunction<ComplexD> ComplexDAcscFunction = new UnaryFunction<ComplexD>(MathEx.Acsc);
		#endregion

		#region Hyperbolic Functions
		public static readonly UnaryFunction<double> DoubleSinhFunction = new UnaryFunction<double>(MathEx.Sinh);
		public static readonly UnaryFunction<double> DoubleCoshFunction = new UnaryFunction<double>(MathEx.Cosh);
		public static readonly UnaryFunction<double> DoubleTanhFunction = new UnaryFunction<double>(MathEx.Tanh);
		public static readonly UnaryFunction<double> DoubleCothFunction = new UnaryFunction<double>(MathEx.Coth);
		public static readonly UnaryFunction<double> DoubleSechFunction = new UnaryFunction<double>(MathEx.Sech);
		public static readonly UnaryFunction<double> DoubleCschFunction = new UnaryFunction<double>(MathEx.Csch);

		public static readonly UnaryFunction<ComplexD> ComplexDSinhFunction = new UnaryFunction<ComplexD>(MathEx.Sinh);
		public static readonly UnaryFunction<ComplexD> ComplexDCoshFunction = new UnaryFunction<ComplexD>(MathEx.Cosh);
		public static readonly UnaryFunction<ComplexD> ComplexDTanhFunction = new UnaryFunction<ComplexD>(MathEx.Tanh);
		public static readonly UnaryFunction<ComplexD> ComplexDCothFunction = new UnaryFunction<ComplexD>(MathEx.Coth);
		public static readonly UnaryFunction<ComplexD> ComplexDSechFunction = new UnaryFunction<ComplexD>(MathEx.Sech);
		public static readonly UnaryFunction<ComplexD> ComplexDCschFunction = new UnaryFunction<ComplexD>(MathEx.Csch);
		#endregion

		#region Hyperbolic Area Functions
		public static readonly UnaryFunction<double> DoubleAsinhFunction = new UnaryFunction<double>(MathEx.Asinh);
		public static readonly UnaryFunction<double> DoubleAcoshFunction = new UnaryFunction<double>(MathEx.Acosh);
		public static readonly UnaryFunction<double> DoubleAtanhFunction = new UnaryFunction<double>(MathEx.Atanh);
		public static readonly UnaryFunction<double> DoubleAcothFunction = new UnaryFunction<double>(MathEx.Acoth);
		public static readonly UnaryFunction<double> DoubleAsechFunction = new UnaryFunction<double>(MathEx.Asech);
		public static readonly UnaryFunction<double> DoubleAcschFunction = new UnaryFunction<double>(MathEx.Acsch);

		public static readonly UnaryFunction<ComplexD> ComplexDAsinhFunction = new UnaryFunction<ComplexD>(MathEx.Asinh);
		public static readonly UnaryFunction<ComplexD> ComplexDAcoshFunction = new UnaryFunction<ComplexD>(MathEx.Acosh);
		public static readonly UnaryFunction<ComplexD> ComplexDAtanhFunction = new UnaryFunction<ComplexD>(MathEx.Atanh);
		public static readonly UnaryFunction<ComplexD> ComplexDAcothFunction = new UnaryFunction<ComplexD>(MathEx.Acoth);
		public static readonly UnaryFunction<ComplexD> ComplexDAsechFunction = new UnaryFunction<ComplexD>(MathEx.Asech);
		public static readonly UnaryFunction<ComplexD> ComplexDAcschFunction = new UnaryFunction<ComplexD>(MathEx.Acsch);
		#endregion

		#region Abs Functions
		/// <summary>
		/// Absolute value function for single-precision floating point numbers.
		/// </summary>
		public static readonly UnaryFunction<float> FloatAbsFunction = new UnaryFunction<float>(System.Math.Abs);
		/// <summary>
		/// Absolute value function for double-precision floating point numbers.
		/// </summary>
		public static readonly UnaryFunction<double> DoubleAbsFunction = new UnaryFunction<double>(System.Math.Abs);
		/// <summary>
		/// Absolute value function for integers.
		/// </summary>
		public static readonly UnaryFunction<int> IntAbsFunction = new UnaryFunction<int>(System.Math.Abs);
		#endregion



		public static readonly UnaryFunction<double> DoubleSqrtFunction = new UnaryFunction<double>(System.Math.Sqrt);
		#endregion

		#region Abs
		/// <summary>
		/// Calculates the absolute value of an integer.
		/// </summary>
		/// <param name="x">An integer.</param>
		/// <returns>The absolute value of <paramref name="x"/>.</returns>
		public static int Abs(int x)
		{
			return System.Math.Abs(x);
		}
		/// <summary>
		/// Calculates the absolute value of a single-precision floating point number.
		/// </summary>
		/// <param name="x">A single-precision floating point number.</param>
		/// <returns>The absolute value of <paramref name="x"/>.</returns>
		public static float Abs(float x)
		{
			return System.Math.Abs(x);
		}
		/// <summary>
		/// Calculates the absolute value of a double-precision floating point number.
		/// </summary>
		/// <param name="x">A double-precision floating point number.</param>
		/// <returns>The absolute value of <paramref name="x"/>.</returns>
		public static double Abs(double x)
		{
			return System.Math.Abs(x);
		}
		/// <summary>
		/// Creates a new array of integers whose element values are the
		/// result of applying the absolute function on the elements of the
		/// given integers array.
		/// </summary>
		/// <param name="array">An array of integers.</param>
		/// <returns>A new array of integers whose values are the result of applying the absolute function to each element in <paramref name="array"/></returns>
		public static int[] Abs(int[] array)
		{
			int[] result = new int[array.Length];
			for (int i = 0; i < array.Length; i++)
			{
				result[i] = Abs(array[i]);
			}
			return result;
		}		
		/// <summary>
		/// Creates a new <see cref="float[]"/> whose element values are the
		/// result of applying the absolute function on the elements of the
		/// given floats array.
		/// </summary>
		/// <param name="array">An array of single-precision floating point values.</param>
		/// <returns>A new <see cref="FloatArrayList"/> whose values are the result of applying the absolute function to each element in <paramref name="array"/></returns>
		public static float[] Abs(float[] array)
		{
			float[] result = new float[array.Length];
			for (int i = 0; i < array.Length; i++)
			{
				result[i] = Abs(array[i]);
			}
			return result;
		}
		/// <summary>
		/// Creates a new <see cref="double[]"/> whose element values are the
		/// result of applying the absolute function on the elements of the
		/// given doubles array.
		/// </summary>
		/// <param name="array">An array of double-precision floating point values.</param>
		/// <returns>A new <see cref="double[]"/> whose values are the result of applying the absolute function to each element in <paramref name="array"/></returns>
		public static double[] Abs(double[] array)
		{
			double[] result = new double[array.Length];
			for (int i = 0; i < array.Length; i++)
			{
				result[i] = Abs(array[i]);
			}

			return result;
		}
		#endregion

		#region Clamp
		/// <summary>
		/// Clamp a <paramref name="value"/> to <paramref name="calmpedValue"/> if it is withon the <paramref name="tolerance"/> range.
		/// </summary>
		/// <param name="value">The value to clamp.</param>
		/// <param name="calmpedValue">The clamped value.</param>
		/// <param name="tolerance">The tolerance value.</param>
		/// <returns>
		/// Returns the clamped value.
		/// result = (tolerance > Abs(value-calmpedValue)) ? calmpedValue : value;
		/// </returns>
		public static double Clamp(double value, double calmpedValue, double tolerance)
		{
			return (tolerance > Abs(value - calmpedValue)) ? calmpedValue : value;
		}
		/// <summary>
		/// Clamp a <paramref name="value"/> to <paramref name="calmpedValue"/> using the default tolerance value.
		/// </summary>
		/// <param name="value">The value to clamp.</param>
		/// <param name="calmpedValue">The clamped value.</param>
		/// <returns>
		/// Returns the clamped value.
		/// result = (EpsilonD > Abs(value-calmpedValue)) ? calmpedValue : value;
		/// </returns>
		/// <remarks><see cref="MathFunctions.EpsilonD"/> is used for tolerance.</remarks>
		public static double Clamp(double value, double calmpedValue)
		{
			return (EpsilonD > Abs(value - calmpedValue)) ? calmpedValue : value;
		}
		/// <summary>
		/// Clamp a <paramref name="value"/> to <paramref name="calmpedValue"/> if it is withon the <paramref name="tolerance"/> range.
		/// </summary>
		/// <param name="value">The value to clamp.</param>
		/// <param name="calmpedValue">The clamped value.</param>
		/// <param name="tolerance">The tolerance value.</param>
		/// <returns>
		/// Returns the clamped value.
		/// result = (tolerance > Abs(value-calmpedValue)) ? calmpedValue : value;
		/// </returns>
		public static float Clamp(float value, float calmpedValue, float tolerance)
		{
			return (tolerance > Abs(value - calmpedValue)) ? calmpedValue : value;
		}

        public static int Clamp(int value, int calmpedValue, float tolerance)
        {
            return (tolerance > Abs(value - calmpedValue)) ? calmpedValue : value;
        }
		/// <summary>
		/// Clamp a <paramref name="value"/> to <paramref name="calmpedValue"/> using the default tolerance value.
		/// </summary>
		/// <param name="value">The value to clamp.</param>
		/// <param name="calmpedValue">The clamped value.</param>
		/// <returns>
		/// Returns the clamped value.
		/// result = (EpsilonF > Abs(value-calmpedValue)) ? calmpedValue : value;
		/// </returns>
		/// <remarks><see cref="MathFunctions.EpsilonF"/> is used for tolerance.</remarks>
		public static float Clamp(float value, float calmpedValue)
		{
			return (EpsilonF > Abs(value - calmpedValue)) ? calmpedValue : value;
		}
        public static int Clamp(int value, int calmpedValue)
        {
            return (EpsilonF > Abs(value - calmpedValue)) ? calmpedValue : value;
        }
		#endregion

		#region Trigonometric Functions
		public static double Sin(double value)
		{
			return System.Math.Sin(value);
		}
		public static float Sin(float value)
		{
			return (float)System.Math.Sin(value);
		}
		public static ComplexD Sin(ComplexD value)
		{
			return ComplexD.Sin(value);
		}

		public static double Cos(double value)
		{
			return System.Math.Cos(value);
		}
		public static float Cos(float value)
		{
			return (float)System.Math.Cos(value);
		}
		public static ComplexD Cos(ComplexD value)
		{
			return ComplexD.Cos(value);
		}

		public static double Tan(double value)
		{
			return System.Math.Tan(value);
		}
		public static float Tan(float value)
		{
			return (float)System.Math.Tan(value);
		}
		public static ComplexD Tan(ComplexD value)
		{
			return ComplexD.Tan(value);
		}

		public static double Cot(double value)
		{
			return 1 / System.Math.Tan(value);
		}
		public static float Cot(float value)
		{
			return (float)(1 / System.Math.Tan(value));
		}
		public static ComplexD Cot(ComplexD value)
		{
			return ComplexD.Cot(value);
		}

		public static double Sec(double value)
		{
			return 1 / System.Math.Cos(value);
		}
		public static float Sec(float value)
		{
			return (float)(1 / System.Math.Cos(value));
		}
		public static ComplexD Sec(ComplexD value)
		{
			return ComplexD.Sec(value);
		}

		public static double Csc(double value)
		{
			return 1 / System.Math.Sin(value);
		}
		public static float Csc(float value)
		{
			return (float)(1 / System.Math.Sin(value));
		}
		public static ComplexD Csc(ComplexD value)
		{
			return ComplexD.Csc(value);
		}
		#endregion

		#region Trigonometric Arcus Functions
		public static double Asin(double value)
		{
			return System.Math.Asin(value);
		}
		public static ComplexD Asin(ComplexD value)
		{
			return ComplexD.Asin(value);
		}

		public static double Acos(double value)
		{
			return System.Math.Acos(value);
		}
		public static ComplexD Acos(ComplexD value)
		{
			return ComplexD.Acos(value);
		}

		public static double Atan(double value)
		{
			return System.Math.Atan(value);
		}
		public static ComplexD Atan(ComplexD value)
		{
			return ComplexD.Atan(value);
		}

		public static double Acot(double value)
		{
			return System.Math.Atan(1 / value);
		}
		public static ComplexD Acot(ComplexD value)
		{
			return ComplexD.Acot(value);
		}

		public static double Asec(double value)
		{
			return System.Math.Acos(1 / value);
		}
		public static ComplexD Asec(ComplexD value)
		{
			return ComplexD.Asec(value);
		}

		public static double Acsc(double value)
		{
			return System.Math.Asin(1 / value);
		}
		public static ComplexD Acsc(ComplexD value)
		{
			return ComplexD.Acsc(value);
		}
		#endregion

		#region Hyperbolic Functions 双曲线函数
		public static double Sinh(double value) 
		{
			return System.Math.Sinh(value);
		}
		public static ComplexD Sinh(ComplexD value)
		{
			return ComplexD.Sinh(value);
		}

		public static double Cosh(double value) 
		{
			return System.Math.Cosh(value);
		}
		public static ComplexD Cosh(ComplexD value)
		{
			return ComplexD.Cosh(value);
		}

		public static double Tanh(double value) 
		{
			return System.Math.Tanh(value);
		}
		public static ComplexD Tanh(ComplexD value)
		{
			return ComplexD.Tanh(value);
		}

		public static double Coth(double value) 
		{
			return 1/System.Math.Tanh(value);
		}
		public static ComplexD Coth(ComplexD value)
		{
			return ComplexD.Coth(value);
		}

		public static double Sech(double value) 
		{
			return 1/Cosh(value);
		}
		public static ComplexD Sech(ComplexD value)
		{
			return ComplexD.Sech(value);
		}

		public static double Csch(double value) 
		{
			return 1/Sinh(value);
		}
		public static ComplexD Csch(ComplexD value)
		{
			return ComplexD.Csch(value);
		}
		#endregion

		#region Hyperbolic Area Functions 反双曲线函数
		public static double Asinh(double value)
		{
			return System.Math.Log(value + System.Math.Sqrt(value * value + 1), System.Math.E);
		}
		public static ComplexD Asinh(ComplexD value)
		{
			return ComplexD.Asinh(value);
		}

		public static double Acosh(double value)
		{
			return System.Math.Log(value + System.Math.Sqrt(value - 1) * System.Math.Sqrt(value + 1), System.Math.E);
		}
		public static ComplexD Acosh(ComplexD value)
		{
			return ComplexD.Acosh(value);
		}

		public static double Atanh(double value)
		{
			return 0.5 * System.Math.Log((1 + value) / (1 - value), System.Math.E);
		}
		public static ComplexD Atanh(ComplexD value)
		{
			return ComplexD.Atanh(value);
		}

		public static double Acoth(double value)
		{
			return 0.5 * System.Math.Log((value + 1) / (value - 1), System.Math.E);
		}
		public static ComplexD Acoth(ComplexD value)
		{
			return ComplexD.Acoth(value);
		}

		public static double Asech(double value)
		{
			return Acosh(1 / value);
		}
		public static ComplexD Asech(ComplexD value)
		{
			return ComplexD.Asech(value);
		}

		public static double Acsch(double value)
		{
			return Asinh(1 / value);
		}
		public static ComplexD Acsch(ComplexD value)
		{
			return ComplexD.Acsch(value);
		}
		#endregion

		#region ApproxEquals
		/// <summary>
		/// Tests whether two single-precision floating point numbers are approximately equal using default tolerance value.
		/// </summary>
		/// <param name="a">A single-precision floating point number.</param>
		/// <param name="b">A single-precision floating point number.</param>
		/// <returns><see langword="true"/> if the two vectors are approximately equal; otherwise, <see langword="false"/>.</returns>
		public static bool ApproxEquals(float a, float b)
		{
			return (System.Math.Abs(a - b) <= EpsilonF);
		}
		/// <summary>
		/// Tests whether two single-precision floating point numbers are approximately equal given a tolerance value.
		/// </summary>
		/// <param name="a">A single-precision floating point number.</param>
		/// <param name="b">A single-precision floating point number.</param>
		/// <param name="tolerance">The tolerance value used to test approximate equality.</param>
		/// <returns><see langword="true"/> if the two vectors are approximately equal; otherwise, <see langword="false"/>.</returns>
		public static bool ApproxEquals(float a, float b, float tolerance)
		{
			return (System.Math.Abs(a - b) <= tolerance);
		}
		/// <summary>
		/// Tests whether two double-precision floating point numbers are approximately equal using default tolerance value.
		/// </summary>
		/// <param name="a">A double-precision floating point number.</param>
		/// <param name="b">A double-precision floating point number.</param>
		/// <returns><see langword="true"/> if the two vectors are approximately equal; otherwise, <see langword="false"/>.</returns>
		public static bool ApproxEquals(double a, double b)
		{
			return (System.Math.Abs(a - b) <= EpsilonD);
		}
		/// <summary>
		/// Tests whether two double-precision floating point numbers are approximately equal given a tolerance value.
		/// </summary>
		/// <param name="a">A double-precision floating point number.</param>
		/// <param name="b">A double-precision floating point number.</param>
		/// <param name="tolerance">The tolerance value used to test approximate equality.</param>
		/// <returns><see langword="true"/> if the two vectors are approximately equal; otherwise, <see langword="false"/>.</returns>
		public static bool ApproxEquals(double a, double b, double tolerance)
		{
			return (System.Math.Abs(a - b) <= tolerance);
		}
		#endregion

        #region Sketch Recognise Math Function

        public static float Top(Vector2F[] points,ref int index)
        {
            float top = float.MinValue;
            index = 0;
            for (int i = 0; i < points.Length; i++)
            {
                if (points[i].Y > top)
                {
                    top = points[i].Y;
                    index = i;
                }
            }
            return top;
        }
        public static float Left(Vector2F[] points,ref int index)
        {
            float left = float.MaxValue;
            index = 0;
            for (int i = 0; i < points.Length; i++)
            {
                if (points[i].X < left)
                {
                    left = points[i].X;
                    index = i;
                }
            }
            return left;
        }
        public static float Right(Vector2F[] points,ref int index)
        {

            float right = float.MinValue;
            index = 0;
            for (int i = 0; i < points.Length; i++)
            {
                if (points[i].X > right)
                {
                    right = points[i].X;
                    index = i;
                }
            }
            return right;
        }
        public static float Bottom(Vector2F[] points,ref int index)
        {
            float bottom = float.MaxValue;
            index = 0;
            for (int i = 0; i < points.Length; i++)
            {
                if (points[i].Y < bottom)
                {
                    bottom = points[i].Y;
                    index = i;
                }
            }
            return bottom;
        }
        public static float SumX(Vector2F[] points)
        {
            float sum = 0;
            foreach (var item in points)
            {
                sum += item.X;
            }
            return sum;
        }
        public static float SumX(Stroke stroke)
        {
            return SumX(stroke.Points.ToArray());
        }
        public static float SumX(List<Vector2F> points)
        {
            return SumX(points.ToArray());
        }
        public static float SumY(Vector2F[] points)
        {
            float sum = 0;
            foreach (var item in points)
            {
                sum += item.Y;
            }
            return sum;
        }
        public static float SumY(Stroke stroke)
        {
            return SumY(stroke.Points.ToArray());
        }
        public static float SumY(List<Vector2F> points)
        {
            return SumY(points.ToArray());
        }
       // public static Vector2F 
        /// <summary>
        /// 求平均值
        /// </summary>
        /// <param name="buffer"></param>
        /// <returns></returns>
        public static double Average(double[] buffer)
        {   //Todo UnitTest
            return Average(buffer, buffer.Length);
        }
        /// <summary>
        /// 前N个数求平均值
        /// </summary>
        /// <param name="buffer"></param>
        /// <param name="n"></param>
        /// <returns></returns>
        public static double Average(double[] buffer, int n)
        {   //Todo UnitTest
            System.Diagnostics.Debug.Assert(n <= buffer.Length);
            
            double sum = 0.0;
            double average = 0.0;

            for (int i = 0; i < n; i++)
                sum = sum + buffer[i];
            average = sum / n;
            return average;
        }

        public static int Average(int[] vals)
        {
            int num = 0;
            foreach (int num2 in vals)
            {
                num += num2;
            }
            return (num / vals.Length);
        }
        /// <summary>
        /// 求前N个的最大 最小值 废弃
        /// </summary>
        /// <param name="p"></param>
        /// <param name="n"></param>
        /// <param name="maxvalue"></param>
        /// <param name="minvalue"></param>
        public static void MaxMin_Of_Array(double[] p, int n, ref double maxValue, ref double minValue)
        {   //Todo UnitTest
            System.Diagnostics.Debug.Assert(n <= p.Length);
            double max = 1.0;
            double min = 100000.0;

            for (int i = 0; i < n ; i++)
            {
                if (p[i] > max)
                    max = p[i];
                if (p[i] < min)
                    min = p[i];
            }
            maxValue = max;
            minValue = min;
        }
        /// <summary>
        /// 求最大 最小值
        /// </summary>
        /// <param name="p"></param>
        /// <param name="maxValue"></param>
        /// <param name="minValue"></param>
        public static void MaxMin(double[] p, ref double maxValue, ref double minValue)
        {   //Todo UnitTest
            double max = 1.0;
            double min = 100000.0;
            foreach (double item in p)
            {
                if (item > max)
                    max = item;
                if (item < min)
                    min = item;
            }
            maxValue = max;
            minValue = min;
        }
        public static void MaxMin(List<double> p,ref double maxValue,ref double minValue)
        {
            double max = 1.0;
            double min = 100000.0;
            foreach (double item in p)
            {
                if (item > max)
                    max = item;
                if (item < min)
                    min = item;
            }
            maxValue = max;
            minValue = min;
        }
        /// <summary>
        /// Get a new Normalize Array based on the old one;
        /// Please Make Sure all the data are  bigger than 0 
        /// and different from each other
        /// </summary>
        /// <param name="buffer"></param>
        /// <returns></returns>
        public static double[] GetNormalizedArray(double[] buffer)
        {   
            double max = 0.0, min = 0.0;
            MaxMin(buffer, ref max, ref min);
            int len = buffer.Length;
            double[] pu = new double[len];
            double dis = max - min;

            System.Diagnostics.Debug.Assert(dis != 0);

            for (int i = 0; i < len; i++)
            {
                pu[i] = (buffer[i] - min) / dis;
            }
            return pu;
        }
        /// <summary>
        /// Normalizing the Array
        /// Please Make Sure all the data are  bigger than 0
        /// and different from each other
        /// </summary>
        /// <param name="buffer"></param>
        public static void NormalizingArray(double[] buffer)
        {
            double max = 0.0, min = 0.0;
            MaxMin(buffer, ref max, ref min);
            int len = buffer.Length;
          
            double dis = max - min;

            System.Diagnostics.Debug.Assert(dis != 0);

            for (int i = 0; i < len; i++)
            {
                buffer[i] = (buffer[i] - min) / dis;
            }
        }
        /// <summary>
        /// Normalizing the array in Reverse way
        /// </summary>
        /// <param name="buffer"></param>
        public static void UnNormalizingArray(double[] buffer)
        {
            double max = 0.0, min = 0.0;
            MaxMin(buffer, ref max, ref min);
            int len = buffer.Length;
            System.Diagnostics.Debug.Assert(max != 0);
            
            for (int i = 0; i < len; i++)
            {
                buffer[i] = (max - buffer[i]) / max;
            }
          
        }
        public static void UnNormalizingArray(List<double> buffer)
        {
            double max = 0.0, min = 0.0;
            MaxMin(buffer, ref max, ref min);
            System.Diagnostics.Debug.Assert(max != 0);

            int len = buffer.Count;
            for (int i = 0; i < len; i++)
            {
                buffer[i] = (max - buffer[i]) / max;
            }
        }

        /// <summary>
        /// 四舍五入绝对值
        /// </summary>
        /// <param name="f"></param>
        /// <returns></returns>
        public static int IntC(double f)
        {   //Todo UnitTest
            if (f > 0)
            {
                if (f - Convert.ToInt32(f) > 0.5)
                    return Convert.ToInt32(f) + 1;
                else
                    return Convert.ToInt32(f);
            }
            else
            {
                if (System.Math.Abs(f - Convert.ToInt32(f)) > 0.5)
                    return Convert.ToInt32(f) - 1;
                else
                    return Convert.ToInt32(f);
            }

        }
        /// <summary>
        /// 四舍五入绝对值
        /// </summary>
        /// <param name="d"></param>
        /// <returns></returns>
        public static int Round(double d)
        {
            int num = (int)d;
            double num2 = (d >= 0.0) ? 0.5 : -0.5;
            if ((d - num) > num2)
            {
                return (num + 1);
            }
            return num;
        }
        /// <summary>
        /// 四舍五入绝对值
        /// </summary>
        /// <param name="f"></param>
        /// <returns></returns>
        public static int Round(float f)
        {
            int num = (int)f;
            float num2 = (f >= 0f) ? 0.5f : -0.5f;
            if ((f - num) > num2)
            {
                return (num + 1);
            }
            return num;
        }

        public static double Square(double v)
        {
            return (v * v);
        }

        public static int Square(int v)
        {
            return (v * v);
        }

        public static double Determinant2x2(double a11, double a12, double a21, double a22)
        {
            return ((a11 * a22) - (a12 * a21));
        }

        /// <summary>
        /// 直线拟合
        /// </summary>
        /// <param name="pt"></param>
        /// <param name="k"></param>
        /// <param name="b"></param>
        /// <returns></returns>
        public static bool FitToLine(Vector2F[] pt,ref double k, ref double b)
        {
            double mX = 0, mY = 0, mXX = 0, mXY = 0;

            if (pt.Length < 2)
            {
                //k = b = 0;
                return false;
            }
            int n = pt.Length;
            foreach (Vector2F item in pt)
            {
                mX += item.X;
                mY += item.Y;
                mXX += item.X * item.X;
                mXY += item.X * item.Y;
            }
            if (mX * mX - mXX * n == 0)
            {
                k = Int32.MaxValue;
                b = 0;
                return false;
            }
            else
            {
                k = (mY * mX - mXY * n) / (mX * mX - mXX * n);
                b = (mY - mX * k) / n;
                return true;
            }
        }
        /// <summary>
        /// 直线拟合+误差
        /// </summary>
        /// <param name="pt"></param>
        /// <param name="err"></param>
        /// <param name="k"></param>
        /// <param name="b"></param>
        /// <returns></returns>
        public static bool FitToLine(Vector2F[] pt,ref double err,ref double k,ref double b)
        {
            double mX = 0, mY = 0, mXX = 0, mXY = 0;

            if (pt.Length < 2)
            {
                //k = b = 0;
                return false;
            }
            int n = pt.Length;
            foreach (Vector2F item in pt)
            {
                mX += item.X;
                mY += item.Y;
                mXX += item.X * item.X;
                mXY += item.X * item.Y;
            }
            if (mX * mX - mXX * n == 0) //垂直或共点
            {
                k = Int32.MaxValue;
                b = 0;
                err = 0;
                return false;
            }
            else
            {
                k = (mY * mX - mXY * n) / (mX * mX - mXX * n);
                b = (mY - mX * k) / n;
                #region fit error
                err = FitLineError(pt, k, b);
                #endregion
            }
            
            return true;
        }

        public static double FitLineError(Vector2F[] pt, double k, double b)
        {
            double err = 0;

            foreach (var item in pt)
            {
                err += Math.Abs(item.Y - (k * item.X + b));
            }
            err = err / pt.Length;

            return err;
        }


        /// <summary>
        /// Bezier Spline methods
        /// </summary>
        public static class BezierSpline
        {

            /// <summary>
            /// Get Bezier Spline Points from Knots;
            /// [knot,control1,control2,knot,control1,control2,knot,control1...knot]
            /// </summary>
            /// <param name="knots"></param>
            /// <returns>Bezier Spline Points</returns>
            public static Vector2F[] GetBezierSplinePoints(Vector2F[] knots)
            {
                Vector2F[] firstControlPoints, secondControlPoints;
                GetCurveControlPoints(knots,out firstControlPoints,out secondControlPoints);
                Vector2F[] bezierSplinePoints = new Vector2F[knots.Length + firstControlPoints.Length + secondControlPoints.Length];
                
                int index = 0;
                bezierSplinePoints[index++] = knots[0];
                for (int i = 0; i < firstControlPoints.Length; ++i)
                {
                    bezierSplinePoints[index++] = firstControlPoints[i];
                    bezierSplinePoints[index++] = secondControlPoints[i];
                    bezierSplinePoints[index++] = knots[i + 1];
                   // lines.Add(new BezierSegment(firstControlPoints[i], secondControlPoints[i], points[i + 1], true));
                }
                return bezierSplinePoints; ;
            }
            /// <summary>
            /// Get open-ended Bezier Spline Control Points.
            /// </summary>
            /// <param name="knots">Input Knot Bezier spline points.</param>
            /// <param name="firstControlPoints">Output First Control points array of knots.Length - 1 length.</param>
            /// <param name="secondControlPoints">Output Second Control points array of knots.Length - 1 length.</param>
            public static void GetCurveControlPoints(Vector2F[] knots, out Vector2F[] firstControlPoints, out Vector2F[] secondControlPoints)
            {
                int n = knots.Length - 1;
                if (n < 1)
                {
                    firstControlPoints = new Vector2F[0];
                    secondControlPoints = new Vector2F[0];
                    return;
                }

                // Calculate first Bezier control points
                // Right hand side vector
                double[] rhs = new double[n];

                // Set right hand side X values
                for (int i = 1; i < n - 1; ++i)
                    rhs[i] = 4 * knots[i].X + 2 * knots[i + 1].X;
                rhs[0] = knots[0].X + 2 * knots[1].X;
                rhs[n - 1] = 3 * knots[n - 1].X;
                // Get first control points X-values
                double[] x = GetFirstControlPoints(rhs);

                // Set right hand side Y values
                for (int i = 1; i < n - 1; ++i)
                    rhs[i] = 4 * knots[i].Y + 2 * knots[i + 1].Y;
                rhs[0] = knots[0].Y + 2 * knots[1].Y;
                rhs[n - 1] = 3 * knots[n - 1].Y;
                // Get first control points Y-values
                double[] y = GetFirstControlPoints(rhs);

                // Fill output arrays.
                firstControlPoints = new Vector2F[n];
                secondControlPoints = new Vector2F[n];
                for (int i = 0; i < n; ++i)
                {
                    // First control point
                    firstControlPoints[i] = new Vector2F(x[i], y[i]);
                    // Second control point
                    if (i < n - 1)
                        secondControlPoints[i] = new Vector2F(2 * knots[i + 1].X - x[i + 1], 2 * knots[i + 1].Y - y[i + 1]);
                    else
                        secondControlPoints[i] = new Vector2F((knots[n].X + x[n - 1]) / 2, (knots[n].Y + y[n - 1]) / 2);
                }
            }

            
            /// <summary>
            /// Get open-ended Bezier Spline Control Points.
            /// </summary>
            /// <param name="knots">Input Knot Bezier spline points.</param>
            /// <param name="firstControlPoints">Output First Control points array of knots.Length - 1 length.</param>
            /// <param name="secondControlPoints">Output Second Control points array of knots.Length - 1 length.</param>
            public static void GetCurveControlPoints(PointF[] knots, out PointF[] firstControlPoints, out PointF[] secondControlPoints)
            {
                int n = knots.Length - 1;
                if (n < 1)
                {
                    firstControlPoints = new PointF[0];
                    secondControlPoints = new PointF[0];
                    return;
                }

                // Calculate first Bezier control points
                // Right hand side vector
                double[] rhs = new double[n];

                // Set right hand side X values
                for (int i = 1; i < n - 1; ++i)
                    rhs[i] = 4 * knots[i].X + 2 * knots[i + 1].X;
                rhs[0] = knots[0].X + 2 * knots[1].X;
                rhs[n - 1] = 3 * knots[n - 1].X;
                // Get first control points X-values
                double[] x = GetFirstControlPoints(rhs);

                // Set right hand side Y values
                for (int i = 1; i < n - 1; ++i)
                    rhs[i] = 4 * knots[i].Y + 2 * knots[i + 1].Y;
                rhs[0] = knots[0].Y + 2 * knots[1].Y;
                rhs[n - 1] = 3 * knots[n - 1].Y;
                // Get first control points Y-values
                double[] y = GetFirstControlPoints(rhs);

                // Fill output arrays.
                firstControlPoints = new PointF[n];
                secondControlPoints = new PointF[n];
                for (int i = 0; i < n; ++i)
                {
                    // First control point
                    firstControlPoints[i] = new PointF((float)x[i], (float)y[i]);
                    // Second control point
                    if (i < n - 1)
                        secondControlPoints[i] = new PointF((float)(2 * knots[i + 1].X - x[i + 1]),(float)(2 * knots[i + 1].Y - y[i + 1]));
                    else
                        secondControlPoints[i] = new PointF((float)((knots[n].X + x[n - 1]) / 2), (float)((knots[n].Y + y[n - 1]) / 2));
                }
            }

            /// <summary>
            /// Solves a tridiagonal system for one of coordinates (x or y) of first Bezier control points.
            /// </summary>
            /// <param name="rhs">Right hand side vector.</param>
            /// <returns>Solution vector.</returns>
            private static double[] GetFirstControlPoints(double[] rhs)
            {
                int n = rhs.Length;
                double[] x = new double[n]; // Solution vector.
                double[] tmp = new double[n]; // Temp workspace.

                double b = 2.0;
                x[0] = rhs[0] / b;
                for (int i = 1; i < n; i++) // Decomposition and forward substitution.
                {
                    tmp[i] = 1 / b;
                    b = (i < n - 1 ? 4.0 : 2.0) - tmp[i];
                    x[i] = (rhs[i] - x[i - 1]) / b;
                }
                for (int i = 1; i < n; i++)
                    x[n - i - 1] -= tmp[n - i] * x[n - i]; // Backsubstitution.

                return x;
            }
        }
       
        #endregion

        #region 弧长 转 角度
        /// <summary>
        /// 单位 弧长 转 角度
        /// </summary>
        /// <param name="radians">单位</param>
        /// <returns></returns>
        public static double RadiansToDegrees(double radians)
        {
            return radians * (180 / Math.PI);
        }
        /// <summary>
        /// 单位角度 转 弧长
        /// </summary>
        /// <param name="degrees"></param>
        /// <returns></returns>
        public static double DegressToRadians(double degrees)
        {
            return degrees * (Math.PI / 180);
        }
        #endregion

        #region 计算几何
        /// <summary>
        /// When All the points are same,there is a bug!
        /// 逻辑坐标 顺时针
        /// </summary>
        /// <param name="pts"></param>
        /// <returns></returns>
        public static Vector2F[] ConvexHull(Vector2F[] pts)
        {
            CPoint[] points = new CPoint[pts.Length];
            for (int i = 0; i < pts.Length; i++)
            {
                points[i] = new CPoint(pts[i]);
            }
            points = convexhull(points);
            Vector2F[] rpoints = new Vector2F[points.Length];
            for (int i = 0; i < rpoints.Length; i++)
            {
                rpoints[i] = points[i].ToVector2F();
            }
            return rpoints;
        }
        public static Vector2F[] ConvexHull(Stroke stroke)
        {
            return ConvexHull(stroke.Points);
        }
        public static Vector2F[] ConvexHull(List<Vector2F> points)
        {
            CPoint[] pts = new CPoint[points.Count];
            for (int i = 0; i < points.Count; i++)
            {
                pts[i] = new CPoint(points[i]);
            }
            pts = convexhull(pts);
            Vector2F[] rpoints = new Vector2F[pts.Length];
            for (int i = 0; i < rpoints.Length; i++)
            {
                rpoints[i] = pts[i].ToVector2F();
            }
            return rpoints;
        }

        private static CPoint[] convexhull(CPoint[] pts)
        {
            // Sort points lexicographically by increasing (x, y)
            int N = pts.Length;
            if (N <= 1)
                return pts;
            Polysort.Quicksort_L(pts);
            CPoint left = pts[0];
            CPoint right = pts[N - 1];
            System.Diagnostics.Debug.WriteLine(left.ToString() + " " + right.ToString());
            // Partition into lower hull and upper hull
            CDLL lower = new CDLL(left), upper = new CDLL(left);
            for (int i = 0; i < N; i++)
            {
                double det = CPoint.Area2(left, right, pts[i]);
                if (det > 0)
                    upper = upper.Append(new CDLL(pts[i]));
                else if (det < 0)
                    lower = lower.Prepend(new CDLL(pts[i]));
            }
            lower = lower.Prepend(new CDLL(right));
            upper = upper.Append(new CDLL(right)).Next;
            // Eliminate points not on the hull
            eliminate(lower);
            eliminate(upper);
            // Eliminate duplicate endpoints
            if (lower.Prev.val.Equals(upper.val))
                lower.Prev.Delete();
            if (upper.Prev.val.Equals(lower.val))
                upper.Prev.Delete();
            // Join the lower and upper hull
            CPoint[] res = new CPoint[lower.Size() + upper.Size()];
            lower.CopyInto(res, 0);
            upper.CopyInto(res, lower.Size());
            return res;

        }

        private static void eliminate(CDLL start)
        {
            CDLL v = start, w = start.Prev;
            bool fwd = false;
            while (v.Next != start || !fwd)
            {
                if (v.Next == w)
                    fwd = true;
                if (CPoint.Area2(v.val, v.Next.val, v.Next.Next.val) < 0) // right turn
                    v = v.Next;
                else
                {                                       // left turn or straight
                    v.Next.Delete();
                    v = v.Prev;
                }
            }
        }

        #region Transform
        /// <summary>
        /// |m11,m12||x| = (m11*x+m12*y) + (m21*x+m22*y)
        /// |m21,m22||y| 
        /// </summary>
        /// <param name="points"></param>
        /// <param name="matrix"></param>
        public static void Transform(Vector2F[] points,Matrix2F matrix)
        {
            float x;
            float y;
            
            for (int i = 0; i < points.Length; i++)
            {
                x = (matrix.M11 * points[i].X) + (matrix.M12 * points[i].Y);
                y = (matrix.M21 * points[i].X) + (matrix.M22 * points[i].Y);
                points[i].X = x;
                points[i].Y = y; 
            }
        }
        /// <summary>
        /// 
        /// </summary>
        /// <param name="points">old buffer</param>
        /// <param name="matrix"></param>
        /// <param name="newPoints">new buffer</param>
        public static void Transform(Vector2F[] oldPoints, Matrix2F matrix, ref Vector2F[] newPoints)
        {
            newPoints = new Vector2F[oldPoints.Length];
            for (int i = 0; i < oldPoints.Length; i++)
            {
                newPoints[i] = Matrix2F.Transform(matrix, oldPoints[i]);
            }
        }
        /// <summary>
        /// Transform Matrix  and Shift Points 
        /// </summary>
        /// <param name="points"></param>
        /// <param name="matrix">Transform Matrix</param>
        /// <param name="original">Shift Point</param>
        public static void Transform(Vector2F[] points,Matrix2F matrix,Vector2F original)
        {

        }
        /// <summary>
        /// Transform Matrix and Shift Points
        /// </summary>
        /// <param name="oldPoints"></param>
        /// <param name="matrix"></param>
        /// <param name="original"></param>
        /// <param name="newPoints"></param>
        public static void Transform(Vector2F[] oldPoints, Matrix2F matrix, Vector2F original,ref Vector2F[] newPoints)
        { 

        }
        /// <summary>
        /// Shift Points then Transform Matrix on old PointsBuffer 
        /// </summary>
        /// <param name="points"></param>
        /// <param name="original"></param>
        /// <param name="matrix"></param>
        public static void Transform(Vector2F[] points, Vector2F original, Matrix2F matrix)
        {
 
        }
        /// <summary>
        /// Shift Points then Transform Matrix to new PointsBuffer
        /// </summary>
        /// <param name="points"></param>
        /// <param name="original"></param>
        /// <param name="matrix"></param>
        /// <param name="newPoints"></param>
        public static void Transform(Vector2F[] oldPoints, Vector2F original, Matrix2F matrix,ref Vector2F[] newPoints)
        {
            newPoints = new Vector2F[oldPoints.Length];
            for (int i = 0; i < oldPoints.Length; i++)
            {
                newPoints[i] = Matrix2F.Transform(matrix, oldPoints[i]-original);
            }

        }
        /// <summary>
        /// rotate some degree surround some point
        /// </summary>
        /// <param name="oldPoints"></param>
        /// <param name="original">绕某点旋转</param>
        /// <param name="matrix">旋转矩阵</param>
        /// <param name="newPoints"></param>
        public static void Rotation(Vector2F[] oldPoints, Vector2F original, Matrix2F matrix, ref Vector2F[] newPoints)
        {
            newPoints = new Vector2F[oldPoints.Length];
            for (int i = 0; i < oldPoints.Length; i++)
            {
                newPoints[i] = Matrix2F.Transform(matrix, oldPoints[i] - original);
                newPoints[i] = newPoints[i] + original;
            }
        }


        #endregion

        #region Min Area Enclosing Rect
        /// <summary>
        /// 最小外包矩形
        /// </summary>
        /// <param name="points">逆时针有序,凸包点,逻辑坐标</param>
        public static double MinAreaEnclosingRectangle(Vector2F[] points)
        {
            Vector2F[] vertex = null;
            double h =0,w =0;
            return MinAreaEnclosingRectangle(points, ref vertex, ref h, ref w);
        }
        /// <summary>
        /// 
        /// </summary>
        /// <param name="points"></param>
        /// <param name="vertex">时序不清,逻辑坐标 逆时,屏幕坐标,顺时,4个顶点</param>
        /// <returns></returns>
        public static double MinAreaEnclosingRectangle(Vector2F[] points,ref Vector2F[] vertex)
        {
            double height = 0, width = 0;
            return MinAreaEnclosingRectangle(points, ref vertex, ref height, ref width);
        }
        /// <summary>
        /// 
        /// </summary>
        /// <param name="points"></param>
        /// <param name="vertex">逆时序</param>
        /// <param name="Height"></param>
        /// <param name="Width"></param>
        /// <returns></returns>
        public static double MinAreaEnclosingRectangle(Vector2F[] points, ref Vector2F[] vertex, ref double Height, ref double Width)
        {
            if (points.Length <= 2)
            {
                vertex = new Vector2F[4];
                vertex[0] = Vector2F.Zero;
                vertex[1] = Vector2F.Zero;
                vertex[2] = Vector2F.Zero;
                vertex[3] = Vector2F.Zero;
                Height = Width = 0;
                return 0;
            }
            Vector2F sVect, eVect, initVect;
            Vector2F[] newPoints = null;
            Vector2F lefttop = new Vector2F(), righttop = new Vector2F(),
                     leftbottom = new Vector2F(), rightbottom = new Vector2F();
            float length, m11, m12, m21, m22;
            float area = float.MaxValue;
            int index = 0;
            #region
            for (int i = 0; i < points.Length; i++)
            {

                sVect = points[i];
                eVect = (i + 1 == points.Length) ? points[0] : points[i + 1];


                initVect = eVect - sVect;

                length = initVect.GetLength();


                m11 = initVect.X / length;
                m12 = initVect.Y / length;
                m21 = -initVect.Y / length;
                m22 = initVect.X / length;

                //System.Console.WriteLine(Math.Acos(m11));
                //System.Console.WriteLine(MathEx.RadiansToDegrees(Math.Acos(m11)));
                Matrix2F matrix = new Matrix2F(m11, m12, m21, m22);

                Rotation(points, points[i], matrix, ref newPoints);

                float top, bottom, left, right;
                int tIndex = 0, bIndex = 0, lIndex = 0, rIndex = 0;

                top = Top(newPoints, ref tIndex);

                bottom = Bottom(newPoints, ref bIndex);
                left = Left(newPoints, ref lIndex);
                right = Right(newPoints, ref rIndex);


                float _area = Math.Abs(top - bottom) * Math.Abs(left - right);

                if (_area < area)
                {
                    area = _area;
                    index = i;
                    lefttop = new Vector2F(left, top);
                    leftbottom = new Vector2F(left, bottom);
                    righttop = new Vector2F(right, top);
                    rightbottom = new Vector2F(right, bottom);
                    Height = Math.Abs(top - bottom);
                    Width = Math.Abs(left - right);
                }
            }
            #endregion

            sVect = points[index];


            eVect = (index + 1 == points.Length) ? points[0] : points[index + 1];

            initVect = eVect - sVect;

            //负角度旋转
            length = initVect.GetLength();
            m11 = initVect.X / length;
            m12 = -initVect.Y / length;
            m21 = initVect.Y / length;
            m22 = initVect.X / length;


            Matrix2F rmatrix = new Matrix2F(m11, m12, m21, m22);

            Vector2F[] nvertex = new Vector2F[4];
            vertex = new Vector2F[4];
            nvertex[0] = rightbottom;
            nvertex[1] = righttop;
            nvertex[2] = lefttop;
            nvertex[3] = leftbottom;
            //Todo 输出old 和角度
            Rotation(nvertex, points[index], rmatrix, ref vertex);
            return area;
        }
      
        public static double MinAreaEnRectangle(List<Vector2F> points)
        {
            return 0;
        }
        public static double MinAreaEnRectangle(Stroke stroke)
        {
            return 0;
        }
        #endregion
        #region Compute Area of Convex Hull Points
        public static double AreaConvexHull(Vector2F[] chpts)
        {
            double area = 0;

            for (int i = 0; i < chpts.Length-1; i++)
            {
                area += chpts[i].X * chpts[i + 1].Y - chpts[i].Y * chpts[i + 1].X;
                
            }
            area += chpts[chpts.Length - 1].X * chpts[0].Y - chpts[chpts.Length - 1].Y * chpts[0].X;
            return Math.Abs(area / 2);
        }
        public static double AreaConvexHull(List<Vector2F> chpts)
        {
            double area = 0;

            for (int i = 0; i < chpts.Count - 1; i++)
            {
                area += chpts[i].X * chpts[i + 1].Y - chpts[i].Y * chpts[i + 1].X;

            }
            area += chpts[chpts.Count - 1].X * chpts[0].Y - chpts[chpts.Count - 1].Y * chpts[0].X;
            return Math.Abs(area / 2);
        }

        public static double AreaTriangle(Vector2F A,Vector2F B,Vector2F C)
        {
            Vector2F[] triangle = new Vector2F[3] { A, B, C };
            return AreaConvexHull(triangle);
        }
        public static double AreaQuadrilateral(Vector2F A, Vector2F B, Vector2F C, Vector2F D)
        {
            Vector2F[] quadrilateral = new Vector2F[4] { A, B, C, D };
            return AreaConvexHull(quadrilateral);
        }
        #endregion

        #region Compute Area of Largest Triangle in ConvexHull
        public static double AltConvexHull(Vector2F[] chpts)
        {
            Vector2F A = Vector2F.Zero, B = Vector2F.Zero, C = Vector2F.Zero;
            return AltConvexHull(chpts,ref A, ref B, ref C);
        }
        /// <summary>
        /// 同在一条直线上
        /// </summary>
        /// <param name="points"></param>
        /// <param name="triangle"></param>
        /// <returns></returns>
        public static double AltConvexHull(Vector2F[] chpts, ref Vector2F[] triangle)
        {
            triangle = new Vector2F[3];
            triangle.Initialize();
            return AltConvexHull(chpts, ref triangle[0], ref triangle[1], ref triangle[2]);
        }
        public static double AltConvexHull(Vector2F[] chpts, ref Vector2F A, ref Vector2F B, ref Vector2F C)
        {
            if (chpts.Length < 3)
            {
                A = Vector2F.Zero; B = Vector2F.Zero; C = Vector2F.Zero;
                return 0;
            }
            Vector2F AB = Vector2F.Zero; Vector2F AC = Vector2F.Zero;
            A = Vector2F.Zero; B = Vector2F.Zero;C = Vector2F.Zero;
            double area = Double.MinValue, _area = 0;
            for (int i = 0; i < chpts.Length - 2; i++)
            {
                for (int j = i + 1; j < chpts.Length - 1; j++)
                {
                    AB.X = chpts[j].X - chpts[i].X; AB.Y = chpts[j].Y - chpts[i].Y;
                    for (int k = j + 1; k < chpts.Length; k++)
                    {
                        AC.X = chpts[k].X - chpts[i].X; AC.Y = chpts[k].Y - chpts[i].Y;
                        _area = Math.Abs(Vector2F.KrossProduct(AC, AB));
                        if (_area > area)
                        {
                            area = _area;
                            A = chpts[i]; B = chpts[j]; C = chpts[k];
                        }
                    }
                }
            }
            return area / 2;
        }
        #endregion


        #region Compute Area of Largest Quadrilateral
        /// <summary>
        /// 有序
        /// </summary>
        /// <param name="chpts"></param>
        /// <param name="A"></param>
        /// <param name="B"></param>
        /// <param name="C"></param>
        /// <param name="D"></param>
        /// <returns></returns>
        public static double AlqConvexHull(Vector2F[] chpts,ref Vector2F A,ref Vector2F B,ref Vector2F C,ref Vector2F D)
        {
            A = Vector2F.Zero; B = Vector2F.Zero; C = Vector2F.Zero; D = Vector2F.Zero;
            if (chpts.Length < 4)
            {
                return 0;
            }
            double area = Double.MinValue, _area = 0;
            for (int i = 0; i < chpts.Length - 3; i++)
            {
                for (int j = i + 1; j < chpts.Length - 2; j++)
                {
                    
                    for (int k = j + 1; k < chpts.Length-1; k++)
                    {
                        for (int t = k+1; t < chpts.Length; t++)
                        {

                            _area = MathEx.Abs(MathEx.AreaQuadrilateral(chpts[i],chpts[j],chpts[k],chpts[t]));
                            if (_area > area)
                            {
                                area = _area;
                                A = chpts[i]; B = chpts[j]; C = chpts[k]; D = chpts[t];
                            }
                        }
                        
                    }
                }
            }
            return area;
        }
        public static double AlqConvexHull(Vector2F[] chpts,ref Vector2F[] lqPoints)
        {
            lqPoints = new Vector2F[4];
            lqPoints.Initialize();
            return AlqConvexHull(chpts, ref lqPoints[0], ref lqPoints[1], ref lqPoints[2], ref lqPoints[3]);
        }
        #endregion
        /// <summary>
        /// 
        /// </summary>
        /// <typeparam name="T"></typeparam>
        /// <param name="lhs"></param>
        /// <param name="rhs"></param>


        #region
        private static bool LinesIntersect(Vector2F lp1, Vector2F lp2, Vector2F lp3, Vector2F lp4, ref double x, ref double y,
            bool returnpoint,
            bool extendA,
            bool extendB)
        {
            // http://local.wasp.uwa.edu.au/~pbourke/geometry/lineline2d/
            // line a is given by P1 and P2, point of intersect for line a (Pa) and b (Pb)
            // Pa = P1 + ua ( P2 - P1 )
            // Pb = P3 + ub ( P4 - P3 )

            // ua(x) = ub(x) and ua(y) = ub (y)
            // x1 + ua (x2 - x1) = x3 + ub (x4 - x3)
            // y1 + ua (y2 - y1) = y3 + ub (y4 - y3)

            // ua = ((x4-x3)(y1-y3) - (y4-y3)(x1-x3)) / ((x4-x3)(x2-x1) - (x4-x3)(y2-y1))
            // ub = ((x2-x1)(y1-y3) - (y2-y1)(x1-x3)) / ((y4-y3)(x2-x1) - (x4-x3)(y2-y1))

            // intersect point x = x1 + ua (x2 - x1)
            // intersect point y = y1 + ua (y2 - y1) 
            double x1 = lp1.X;
            double x2 = lp2.X;
            double x3 = lp3.X;
            double x4 = lp4.X;
            double y1 = lp1.Y;
            double y2 = lp2.Y;
            double y3 = lp3.Y;
            double y4 = lp4.Y;

            double denominator = ((y4 - y3) * (x2 - x1) - (x4 - x3) * (y2 - y1));
            if (denominator == 0) // lines are parallel
                return false;
            double numerator_ua = ((x4 - x3) * (y1 - y3) - (y4 - y3) * (x1 - x3));
            double numerator_ub = ((x2 - x1) * (y1 - y3) - (y2 - y1) * (x1 - x3));
            double ua = numerator_ua / denominator;
            double ub = numerator_ub / denominator;
            // if a line is not extended then ua (or ub) must be between 0 and 1
            if (extendA == false)
            {
                if (ua < 0 || ua > 1)
                    return false;
            }
            if (extendB == false)
            {
                if (ub < 0 || ub > 1)
                    return false;
            }
            if (extendA || extendB) // no need to chck range of ua and ub if check is one on lines 
            {
                x = x1 + ua * (x2 - x1);
                y = y1 + ua * (y2 - y1);
                return true;
            }
            if (ua >= 0 && ua <= 1 && ub >= 0 && ub <= 1)
            {
                if (returnpoint)
                {
                    x = x1 + ua * (x2 - x1);
                    y = y1 + ua * (y2 - y1);
                }
                return true;
            }
            return false;
        }
        public static bool LinesIntersect(Vector2F lp1, Vector2F lp2, Vector2F lp3, Vector2F lp4)
        {
            double x = 0;
            double y = 0;
            return LinesIntersect(lp1, lp2, lp3, lp4, ref x, ref y, false, false, false);
        }
        public static bool LinesIntersect(Segment line1,Segment line2)
        {
            double x = 0;
            double y = 0;
            return LinesIntersect(line1.P1, line1.P0, line2.P1, line2.P0, ref x, ref y, false, false, false);
        }
        public static Vector2F LinesIntersectPoint(Vector2F lp1, Vector2F lp2, Vector2F lp3, Vector2F lp4)
        {
            double x = 0;
            double y = 0;
            if (LinesIntersect(lp1, lp2, lp3, lp4, ref x, ref y, true, false, false))
                return new Vector2F(x, y);
            return Vector2F.Zero;
        }
        #endregion

        #endregion
        public static void Swap<T>(ref T lhs, ref T rhs)
		{
			T temp = lhs;
			lhs = rhs;
			rhs = temp;
		}

        public static Vector2F[] Point2DSort_Clockwise(Vector2F[] points)
        {
            CPoint[] pts = new CPoint[points.Length];
            for (int i = 0; i < points.Length; i++)
            {
                pts[i] = new CPoint(points[i]);
            }
            Polysort.Quicksort_L(pts);
            int N = pts.Length;
            CPoint left = pts[0];
            CPoint right = pts[N - 1];
            
            // Partition into lower hull and upper hull
            //CDLL lower = new CDLL(left), upper = new CDLL(left);
            List<CPoint> lower = new List<CPoint>();
            List<CPoint> upper = new List<CPoint>();

            lower.Add(left);
            
            upper.Add(right);
           

            for (int i = 0; i < N; i++)
            {
                double det = CPoint.Area2(left, right, pts[i]);
                if (det > 0)
                    upper.Add(pts[i]);
                else if (det < 0)
                    lower.Add(pts[i]);
            }


            CPoint[] l = lower.ToArray();
            CPoint[] u = upper.ToArray();

            Polysort.Quicksort_L(l);
            Polysort.Quicksort_B(u);

            CPoint[] r = new CPoint[l.Length + u.Length];
            l.CopyTo(r, 0);
            u.CopyTo(r, l.Length);

            Vector2F[] repoints = new Vector2F[r.Length];
            for (int i = 0; i < r.Length; i++)
            {
                repoints[i] = new Vector2F(r[i].x, r[i].y);
            }
            return repoints;
        }
        public static Vector2F[] Point2DSort_AntiClockwise(Vector2F[] points)
        {
            CPoint[] pts = new CPoint[points.Length];
            for (int i = 0; i < points.Length; i++)
            {
                pts[i] = new CPoint(points[i]);
            }
            Polysort.Quicksort_B(pts);
            int N = pts.Length;
            CPoint left = pts[0];
            CPoint right = pts[N - 1];

            // Partition into lower hull and upper hull
            //CDLL lower = new CDLL(left), upper = new CDLL(left);
            List<CPoint> lower = new List<CPoint>();
            List<CPoint> upper = new List<CPoint>();

            lower.Add(left);

            upper.Add(right);


            for (int i = 0; i < N; i++)
            {
                double det = CPoint.Area2(left, right, pts[i]);
                if (det > 0)
                    upper.Add(pts[i]);
                else if (det < 0)
                    lower.Add(pts[i]);
            }


            CPoint[] l = lower.ToArray();
            CPoint[] u = upper.ToArray();

            Polysort.Quicksort_B(l);
            Polysort.Quicksort_L(u);

            CPoint[] r = new CPoint[l.Length + u.Length];
            l.CopyTo(r, 0);
            u.CopyTo(r, l.Length);

            Vector2F[] repoints = new Vector2F[r.Length];
            for (int i = 0; i < r.Length; i++)
            {
                repoints[i] = new Vector2F(r[i].x, r[i].y);
            }
            return repoints;
        }
    }


    
    #region AUX CLASS
    class CPoint : Ordered
    {
        private static readonly Random rnd = new Random();

        public double x, y;

        public int X
        {
            get
            {
                return (int)x;
            }
            set
            {
                x = (double)value;
            }
        }
        public int Y
        {
            get
            {
                return (int)y;
            }
            set
            {
                y = (double)value;
            }
        }

        public CPoint(double x, double y)
        {
            this.x = x; this.y = y;
        }
        public CPoint(Vector2F pt)
        {
            this.x = pt.X; this.y = pt.Y;
        }
        public override string ToString()
        {
            return "(" + x + ", " + y + ")";
        }

        public static CPoint Random(int w, int h)
        {
            return new CPoint(rnd.Next(w), rnd.Next(h));
        }
        public Vector2F ToVector2F()
        {
            return new Vector2F(this.X, this.Y);
        }

        public bool Equals(CPoint p2)
        {
            return x == p2.x && y == p2.y;
        }

        public override bool Less(Ordered o2)
        {
            CPoint p2 = (CPoint)o2;
            return x < p2.x || x == p2.x && y < p2.y;
        }
        public override bool Bigger(Ordered o2)
        {
            CPoint p2 = (CPoint)o2;
            return x > p2.x || x == p2.x && y > p2.y;
        }
       

        // Twice the signed area of the triangle (p0, p1, p2)
        public static double Area2(CPoint p0, CPoint p1, CPoint p2)
        {
            return p0.x * (p1.y - p2.y) + p1.x * (p2.y - p0.y) + p2.x * (p0.y - p1.y);
        }
       
    }


    class CDLL
    {
        private CDLL prev, next;     // not null, except in deleted elements
        public CPoint val;

        // A new CDLL node is a one-element circular list
        public CDLL(CPoint val)
        {
            this.val = val; next = prev = this;
        }

        public CDLL Prev
        {
            get { return prev; }
        }

        public CDLL Next
        {
            get { return next; }
        }

        // Delete: adjust the remaining elements, make this one point nowhere
        public void Delete()
        {
            next.prev = prev; prev.next = next;
            next = prev = null;
        }

        public CDLL Prepend(CDLL elt)
        {
            elt.next = this; elt.prev = prev; prev.next = elt; prev = elt;
            return elt;
        }

        public CDLL Append(CDLL elt)
        {
            elt.prev = this; elt.next = next; next.prev = elt; next = elt;
            return elt;
        }

        public int Size()
        {
            int count = 0;
            CDLL node = this;
            do
            {
                count++;
                node = node.next;
            } while (node != this);
            return count;
        }

        public void PrintFwd()
        {
            CDLL node = this;
            do
            {
                Console.WriteLine(node.val);
                node = node.next;
            } while (node != this);
            Console.WriteLine();
        }

        public void CopyInto(CPoint[] vals, int i)
        {
            CDLL node = this;
            do
            {
                vals[i++] = node.val;	// still, implicit checkcasts at runtime 
                node = node.next;
            } while (node != this);
        }
       
    }

    class Polysort
    {
        private static void swap(CPoint[] arr, int s, int t)
        {
            CPoint tmp = arr[s]; arr[s] = arr[t]; arr[t] = tmp;
        }
        private static void swap(Ordered[] arr, int s, int t)
        {
            Ordered tmp = arr[s]; arr[s] = arr[t]; arr[t] = tmp;
        }


        // Typed OO-style quicksort a la Hoare/Wirth

        private static void qsort_l(Ordered[] arr, int a, int b)
        {
            // sort arr[a..b]
            if (a < b)
            {
                int i = a, j = b;
                CPoint x = (CPoint)(arr[(i + j) / 2]);
                do
                {
                    while (arr[i].Less(x)) i++;
                    while (x.Less(arr[j])) j--;
                    if (i <= j)
                    {
                        swap(arr, i, j);
                        i++; j--;
                    }
                } while (i <= j);
                qsort_l(arr, a, j);
                qsort_l(arr, i, b);
            }
        }
        private static void qsort_b(Ordered[] arr, int a, int b)
        {
            if (a < b)
            {
                int i = a, j = b;
                CPoint x = (CPoint)(arr[(i + j) / 2]);
                do
                {
                    while (arr[i].Bigger(x)) i++;
                    while (x.Bigger(arr[j])) j--;
                    if (i <= j)
                    {
                        swap(arr, i, j);
                        i++; j--;
                    }
                } while (i <= j);
                qsort_b(arr, a, j);
                qsort_b(arr, i, b);
            }
        }

        public static void Quicksort_L(Ordered[] arr)
        {
            qsort_l(arr, 0, arr.Length - 1);
        }
        public static void Quicksort_B(Ordered[] arr)
        {
            qsort_b(arr, 0, arr.Length - 1);
        }
    }

    abstract class Ordered
    {
        public abstract bool Less(Ordered that);
        public abstract bool Bigger(Ordered that);
    }
    #endregion

    //#region NUnit Test
    ///// <summary>
    ///// Test MathEx Functions
    ///// </summary>
    //[TestFixture]
    //public class Test_MathEx_Class
    //{
    //    /// <summary>
    //    /// Test FitLine Function
    //    /// </summary>
    //    /// <remarks>
    //    /// - Desired Result:	Test
    //    /// - Coordinator:      FitLine
    //    /// - Conditions:       Function
    //    /// </remarks>
    //    [Test]
    //    public void Test_FitLine_Function()
    //    {
    //        Vector2F[] pt = new Vector2F[7];
    //        pt[0] = new Vector2F(544, 117);
    //        pt[1] = new Vector2F(545, 117);
    //        pt[2] = new Vector2F(540, 115);
    //        pt[3] = new Vector2F(541, 116);
    //        pt[4] = new Vector2F(543, 117);
    //        pt[5] = new Vector2F(535, 114);
    //        pt[6] = new Vector2F(538, 115);



    //        double err = 0, k = 0, b = 0;
    //        MathEx.FitToLine(pt, ref err,ref k, ref b);
    //        System.Console.Error.WriteLine(err + " " + k + " " + b);
    //        Assert.IsTrue(k > 0);
    //    }
    //    /// <summary>
    //    /// Test Normalize_Array
    //    /// </summary>
    //    /// <remarks>
    //    /// - Desired Result:	Test
    //    /// - Coordinator:      Normalize
    //    /// - Conditions:       Array
    //    /// </remarks>
    //    [Test]
    //    public void Test_Normalize_Array()
    //    {
    //        double[] array = new double[] { -1, -2, -3, 4, 5, 6, 7, 8, 9 };
    //        double[] normalize_array = MathEx.GetNormalizedArray(array);
    //        foreach (var item in normalize_array)
    //        {
    //            System.Console.WriteLine(item.ToString()); 
    //        }
            
    //    }
    //    /// <summary>
    //    /// Test Normalizing Array Function
    //    /// </summary>
    //    /// <remarks>
    //    /// - Desired Result:	Test
    //    /// - Coordinator:      NormalizeingArray
    //    /// - Conditions:       Function
    //    /// </remarks>
    //    [Test]
    //    public void Test_NormalizeingArray_Function()
    //    {
    //        double[] array = new double[] { 1, 2, 3, 4, 5, 6, 7, 8, 9 };
    //        MathEx.NormalizingArray(array);
    //        foreach (var item in array)
    //        {
    //            System.Console.WriteLine(item.ToString()); 
    //        } 
    //    }
    //    /// <summary>
    //    /// Test UnNormalizing Array Function
    //    /// </summary>
    //    /// <remarks>
    //    /// - Desired Result:	Test
    //    /// - Coordinator:      UnNormalizing
    //    /// - Conditions:       Function
    //    /// </remarks>
    //    [Test]
    //    public void Test_UnNormalizing_Function()
    //    {
    //        double[] array = new double[] {  };
    //        MathEx.UnNormalizingArray(array);
    //        foreach (var item in array)
    //        {
    //            System.Console.WriteLine(item.ToString());
    //        } 
    //    }
    //    /// <summary>
    //    /// Test Convex Hull Function
    //    /// </summary>
    //    /// <remarks>
    //    /// - Desired Result:	Test
    //    /// - Coordinator:      ConveHull
    //    /// - Conditions:       Function
    //    /// </remarks>
    //    [Test]
    //    public void Test_ConveHull_Function()
    //    {
    //        Vector2F[] pts = new Vector2F[4]{ new Vector2F(1,1),new Vector2F(1,1),new Vector2F(1,1),new Vector2F(1,1)};
    //        pts = MathEx.ConvexHull(pts);
    //        foreach (var item in pts)
    //        {
    //            System.Console.WriteLine(item.ToString());
    //        }
    //    }
    //    /// <summary>
    //    /// Test AreaCovexHull when only 2 points
    //    /// </summary>
    //    /// <remarks>
    //    /// - Desired Result:	Test
    //    /// - Coordinator:      AreaConvexHull
    //    /// - Conditions:       Function        
    //    /// </remarks>
    //    [Test]
    //    public void Test_AreaConvexHull_Function()
    //    {
    //        Vector2F[] pts = new Vector2F[3] { new Vector2F(1, 1), new Vector2F(4, 4) ,new Vector2F(5,5)};
    //        double area = MathEx.AreaConvexHull(pts);
    //        Console.WriteLine(area);
    //    }

    //    /// <summary>
    //    /// Test Largest Triangle
    //    /// </summary>
    //    /// <remarks>
    //    /// - Desired Result:	Test
    //    /// - Coordinator:      AltConvexHull
    //    /// - Conditions:       Function
    //    /// </remarks>
    //    [Test]
    //    public void Test_AltConvexHull_Function()
    //    {
    //        Vector2F[] pts = new Vector2F[3] { new Vector2F(1, 1), new Vector2F(4, 4), new Vector2F(5, 5) };
    //        Vector2F[] vertex = null;
    //        double area = MathEx.AltConvexHull(pts,ref vertex);
    //        foreach (var item in vertex)
    //        {
    //            Console.WriteLine(item.ToString());
    //            Console.WriteLine(area);
    //        }
    //    }
    //    /// <summary>
    //    /// Test AlqConvexHull Function
    //    /// </summary>
    //    /// <remarks>
    //    /// - Desired Result:	Test
    //    /// - Coordinator:      AlqConvexHull
    //    /// - Conditions:       Function
    //    /// </remarks>
    //    [Test]
    //    public void Test_AlqConvexHull_Function()
    //    {
    //        Vector2F[] pts = new Vector2F[3] { new Vector2F(1, 1), new Vector2F(4, 4), new Vector2F(5, 5) };
    //        Vector2F a = Vector2F.Zero, b = Vector2F.Zero, c = Vector2F.Zero, d = Vector2F.Zero;
    //        double area = MathEx.AlqConvexHull(pts,ref a,ref b,ref c,ref d);
    //        Console.WriteLine(area);
    //    }

    //    /// <summary>
    //    /// Test Line Intersect
    //    /// </summary>
    //    /// <remarks>
    //    /// - Desired Result:	Test
    //    /// - Coordinator:      Line
    //    /// - Conditions:       Intersect
    //    /// </remarks>
    //    [Test]
    //    public void Test_Line_Intersect()
    //    {
    //        Vector2F p1 = new Vector2F(0, 0);
    //        Vector2F p2 = new Vector2F(0, 2);
    //        Vector2F p3 = new Vector2F(0, 1);
    //        Vector2F p4 = new Vector2F(3, 3);

    //        Assert.IsTrue(MathEx.LinesIntersect(p1,p2,p3,p4));
    //    }
    //    enum MyEnum
    //    {
    //        a,
    //        b,
    //        c
    //    }
    //    /// <summary>
    //    /// Write your test story here.
    //    /// </summary>
    //    /// <remarks>
    //    /// - Desired Result:	Test
    //    /// - Coordinator:      GetNormalizedArray
    //    /// - Conditions:       Function
    //    /// </remarks>
    //    [Test]
    //    public void Test_GetNormalizedArray_Function()
    //    {
    //        Console.WriteLine(MathEx.GetNormalizedArray(new double[] { 1, 2, 3, 4, 5 }));
    //    }

    //    /// <summary>
    //    /// 是否为按顺时针 或 逆时针 排序
    //    /// </summary>
    //    [Test]
    //    public void Test_Polysort_Function()
    //    {
    //        CPoint p1 = new CPoint(0, 0);
    //        CPoint p2 = new CPoint(0, 0);
    //        CPoint p3 = new CPoint(0, 0);
    //        CPoint p4 = new CPoint(0, 0);
    //        CPoint p5 = new CPoint(1, 0);
    //        CPoint p6 = new CPoint(0, 1);
    //        CPoint p7 = new CPoint(0.3, 0);
    //        CPoint p8 = new CPoint(1, 0.5);

    //        CPoint[] pts = new CPoint[8] { p1, p2, p3, p4, p5, p6, p7, p8 };
    //        Polysort.Quicksort_L(pts);

    //        foreach (var item in pts)
    //        {
    //            Console.WriteLine(String.Format("x:{0},y:{1}", item.x, item.y));
    //        }
    //        int N = pts.Length;
    //        CPoint left = pts[0];
    //        CPoint right = pts[N - 1];
    //        System.Diagnostics.Debug.WriteLine(left.ToString() + " " + right.ToString());
    //        // Partition into lower hull and upper hull
    //        //CDLL lower = new CDLL(left), upper = new CDLL(left);
    //        List<CPoint> lower = new List<CPoint>();
    //        List<CPoint> upper = new List<CPoint>();

    //        lower.Add(left);
    //        Console.WriteLine(left.x+" "+left.y);
    //        upper.Add(right);
    //        Console.WriteLine(right.x+" " + right.y);

    //        for (int i = 0; i < N; i++)
    //        {
    //            double det = CPoint.Area2(left, right, pts[i]);
    //            if (det > 0)
    //                upper.Add(pts[i]);
    //            else if (det < 0)
    //                lower.Add(pts[i]);
    //        }
           

    //        CPoint[] l = lower.ToArray();
    //        CPoint[] u = upper.ToArray();

    //        Polysort.Quicksort_L(l);
    //        Polysort.Quicksort_B(u);

    //        CPoint[] r = new CPoint[l.Length + u.Length];
    //        l.CopyTo(r, 0);
    //        u.CopyTo(r, l.Length);


    //        foreach (var item in l)
    //        {
    //            Console.WriteLine(String.Format("x:{0},y:{1}", item.x, item.y));
    //        }
    //        Console.WriteLine("-----------------");
    //        foreach (var item in u)
    //        {
    //            Console.WriteLine(String.Format("x:{0},y:{1}", item.x, item.y));
    //        }
    //        Console.WriteLine("----------------------");
    //        foreach (var item in r)
    //        {
    //             Console.WriteLine(String.Format("x:{0},y:{1}", item.x, item.y));
    //        }


            
    //        //upper.PrintFwd();
    //        //lower.PrintFwd();
            
    //        //foreach (var item in lower_t)
    //        //{
    //        //    Console.WriteLine(String.Format("x:{0},y:{1}", item.x, item.y));
    //        //}
    //        //foreach (var item in upper_t)
    //        //{
    //        //    Console.WriteLine(String.Format("x:{0},y:{1}", item.x, item.y));
    //        //}
    //    }
    //}
    

    //#endregion
}
